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ABSTRACT 

We make use of hydrodynamical simulations of the intergalactic medium (IGM) to create model 
quasar absorption spectra. We compare these model spectra with the observed Keck spectra of three 
z > 6.2 quasars with full Gunn-Peterson troughs: SDSS J1148+5251 (z = 6.42), SDSS J1030+0524 
{z = 6.28), and SDSS J1623+3112 (z = 6.22). We fit the probability density distributions (PDFs) 
of the observed Lyman a optical depths (r Q ) with those generated from the simulation, by adopting 
a template for the quasar's intrinsic spectral shape, and exploring a range of values for the size of 
the quasar's surrounding HII region, the volume-weighted mean neutral hydrogen fraction in the 
ambient IGM, xhi, and the quasar's ionizing photon emissivity, Nq. In order to avoid averaging over 
possibly large sightline-to-sightlinc fluctuations in IGM properties, we analyze each observed quasar 
independently. We find the following results for J1148+5251, J1030+0524, and J1623+3112: The 
best -fit sizes Rs are 40, 41, and 29 (comoving) Mpc, respectively. For the later two quasars, the 
value is significantly larger than the radius corresponding to the wavelength at which the quasar's flux 
vanishes. These constraints are tight, with only ~ 10% uncertainties, comparable to those caused by 
redshift-determination errors. The best-fit values of Nq are 2.1, 1.3, and 0.9 x 10 57 s _1 , respectively, 
with a factor of ~ 2 uncertainty in each case. These values are a factor of ~ 5 lower than expected 
from the template spectrum of Telfer et al. (2002). Finally, the best-fit values of Shi are 0.16, 1.0, and 
1.0, respectively. The uncertainty in the case of J1148+5251 is large, and xm is not well constrained. 
However, for both J1030+0524 and J1623+3112 we find a significant lower limit of xm > 0.033. Our 
method is different from previous analyses of the GP absorption spectra of these quasars, and our 
results strengthen the evidence that the rapid end-stage of reionization is occurring near z ~ 6. 

Subject headings: cosmology: theory - early Universe - galaxies: formation - high-redshift - evolution 
- quasars: spectrum 



1. INTRODUCTION 

The epoch of reionization, when the radiation from 
early generations of astrophysical objects ionized the in- 
tergalactic medium (IGM), offers a wealth of information 
about cosmological structure formation and physical pro- 
cesses in the early universe. Only recently have we begun 
to gather clues concerning this epoch. Th e Thomson 
scattering optical depth, r e = 0.09 ± 0.03 (|Page et al.l 
120061: ISpergel et al.ll2006bh . measured from the polariza- 
tion anisotropics in the cosmic microwave background 
(CMB) by the Wilkenson Microwave Anisotropy Probe 
(WMAP), suggests that cosmic reionization began at 
redshift z > 10. However, the current uncertainty of 
this value, and the fact that the optical depth measures 
only an integrated electron column density, means that 
many possible competing reionization scenarios can still 
be accommodated. 

In order to chart the progress of reionization, several 
attempts have been made to measure the neutral frac- 
tion of the IGM at high redshifts (z > 6). The detec- 
tion of dark spectral regions, so-called Gunn-Peterson 
(GP) troughs, in the Lyman line absorption spectra of 
high-redshift quasars discovered in the Sloan Digital Sky 
Survey (SDSS) offers a probe of the late stages of reion- 
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ization. The Lya optical depths associated with such GP 
troughs, and more importantly their accelerated evolu- 
tion at z > 6, suggests that the reionization epoch is end- 
ing at z ~ 6, with a lower limit on the volume-weighted 
IGM neutral fraction, xm(z ~ 6) > 10~ 3 , along at least 
the line of sight ( LOS) to two quasars (|White et al.l l2003: 
iFan et aT]|2006af ). The sharp decline in flux at the edge 
of the GP trough of quasar J1030+0524 can be used to 
infer the stronger constraint of xm{z ~ 6) > 0.2, along 
that LOS (iMesinger fc Haimanll2004h . The observed size 
of the quasars' surrounding HII regions can be used to 
infer a similar, independent lowe r limit (|Wvithe fc Loebl 
12004 IMesinger fc Haiman1l2004D . although this requires 
additional assumptions concerning the quasar and its en- 
vironment, and the inferred neutral fraction scales di- 
rectly with these a ssumptions (jBolton fc HaehnertJl2006t 
iMaselli et al.l 12006). The rapid evolution in the sizes of 
the HII regions for different quasars between z ~ 5.7 
and z ~ 6.4 suggests that the average neutral frac- 
tion increased by a factor of ~ 10 during this interval 
(|Fan et al.ll2006a( ). although there ar e alternative ways to 
accou nt for the steep size-evolution (jBolton fc HaehnaPTJ 
2006). The distribution of spectral sizes of dark gaps 
should help distingui s h these scenario s in the future 
(jGallerani et alj l2006t IFan et all l2006al) . Upper limits 
on the neutral fraction at z ~ 6 have also been suggested 
by other means. The lack of significant evolution in the 
luminosity function (LF) of Lyman a emitters between 
z w 5.7 and z ~ 6.5 was used to infer xm(z ~ 6) < 0.3 
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(|Malhotra fe Rhoadsll200l lHaiman & Cenll2005lh if the 
galaxies are assumed to lie in isolation, or xm{z ~ 6) < 
0.5 (jFurlanetto et al.ll2006l ). if an analytic prescription of 
the clustering of HII bubbles (jFurlanetto et al.ll2004bl ) is 
taken into account. However, a recent, more accurate 
measurement of the LF at z « 6.5 reve aled evidence of a 
decre ase, weakening these upper limits (|Kashikawa et al.l 
2006). Finally, the lack of noticeable absorption due 
to the Lya damping wing in the afterglow spectrum of 
the gamma-ray burst 050904 at z = 6.3 can constr ain 
x m (z ~ 6) < 0.2 along that LOS (|Totani et al.ll2006[ ). 

The preponderance of constraints mentioned above 
suggests that the combination of all existing data is con- 
sistent with Xm{z ~ 6) ~ 0.1. While this is plausible, 
the constraints listed above all rely on various model as- 
sumptions. For example, the conversion of a n observed 
average flux decrement to a neutral fractio n (|Fan et al.l 
[200l IWhite"et~aLl l200l iFan et alJ l2006ah is difficult, 
and requires precise ab-initio knowl edge of the IGM 
densi t y distribution at high-redshifts (|Oh fc Furlanettol 
120051: [Fan et all l2006a[) . All of the other constraints 
mentioned above rely on the IGM density distribution, 
as well, albeit less s ensitively. Techniques involving the 
sizes of HII regions (IWvithe fc Loebll200l iMaselli et all 
2006; IFan et al.ll2006aft are subject to de generacies be- 
tween quasar lifetime, ionizing flux, and the IGM neutral 
fraction. Further complicating matters for the HII bub- 
ble size techniques is the fact that the onset of the GP 
trough need not correspond to the edge of the HII region. 
The GP through just corresponds to the region where 
the flux drops below the detection threshold, and thus 
merel y provides a lower limit on the size of the HI I re- 
gion ^^^^^^^^^^j^l iMasdiiitan MM, see 
also iBolton fc Haehneltl [20061 for a recent, comprehen- 
sive review of these issues). Furthermore, determining 
the number density and evolution of Lya emitting galax- 
ies is com plicated and dep ends on knowledge of galaxy 
clust ering (|Haiman fc Cenll2005t iFurlanetto et al.ll2004ai . 
2006) , abs orber clustering, as well as the large-scale den- 
sity ( e.g. Barkana fc Loebl 2004b, a) and velocity fields 
(e.g. ICen et all |2005| ). which are all unconstrained at 
the relevant high redshifts. 

More generally, the drawback of many of these meth- 
ods is that they rely on the combination of data from 
multiple LOSs to obtain statistical significance. How- 
ever, assuming constancy in the IGM properties along 
disparate LOSs is dangerous at high-redshifts, especially 
approaching the end of reionization, when large sightline- 
to-sightline fluctuations are expected in the density field 
and the ionizing backg rou nd (see, e.g., rece nt reviews by 
iDjorgovski et al.ll2006l and lLidz et al.ll2006h . Indeed, the 
current sample of high-redshift quasars already reveals 
large fl uctuations in the IGM properties along different 
LOSs ([Fan et all l"2006af ) . This would suggest that the 
most rewarding technique would be able to analyze and 
obtain constraints from each source independently from 
other sources. This is statistically less powerful (e.g. 
iMesinger et alJl2~004f ). but the rewards can be rich, po- 
tentially providing additional constraints on the proper- 
ties of individual quasars and their environments, such as 
the quasar's emissivity of ionizing photons, size of its sur- 
rounding HII region, sharpness of the HII region bound- 
ary, the quasars' X-ray emission, as well as a measure of 



the sightline-to-sightline flu ctuations in these quantities 
(|Mesinger fc Haimanll200l . 

The purpose of this paper is to derive these param- 
eters, i.e. the hydrogen neutral fraction, the size of 
surrounding HII region, and the ionizing emissivity, by 
separately modeling the spectra of each of the z > 6.2 
quasars. Specifically, we analyze J1148+5251 (z = 6.42), 
J1030+0524 (z = 6.28) and J1623+3112 (z = 6.22). At 
present, a 4th quasar is known with a full GP trough. 
This source, however, is a broad absorption line (BAL) 
quasar J1048+4637 (z = 6.2), with complex spectral fea- 
tures (Fan, X., private communication) that preclude a 
simple interpretation in terms of IGM properties. We 
have therefore omitted this quasar from our study. 

Our analysis dif f ers f rom our previous work 
(jMesinger fc Haimanl [2004) where we modeled the 
gross transmission features in Lya and Ly/3 close to 
the edge of the HII region surrounding J1030+0524. In 
the present study, we perform a much more detailed 
analysis, by fitting the quasar's intrinsic emission, and 
modeling the distribution of optical depths, pixel-by- 
pixel, blueward of th e Lya line center, fo llowing the 
procedure proposed i n IMesinger et~aU |2004). Unlike in 
IMesinger fc Haimanl (|2004h . here we do not make use 
of the Ly/3 absorption spectrum, except to set a lower 
limit in the allowed size of the HII region surrounding 
J1030+0524, the only one among the three quasars 
that has a noticeable offset between the Lya and Ly/J 
GP troughs. Our reason for ignoring the Ly/3 region is 
that the majority of the pixels we analyze (see Figure [2] 
below) have flux detected in the Lya absorption region, 
which can directly be converted to a Lya optical depth. 
This avoids the need for a statistical modeling of the 
foreground Lya absorption, necessary to convert Ly/3 to 
a Lya optical depth. 

This paper is organized as follows. In §[2j we describe 
our analysis technique, including generating the simu- 
lated absorption (§ 12. ip . fitting for the quasars' intrinsic 
emission spectrum (§ 12. 2| . and the statistical method 
used to compare the observed and simulated spectra 
(§ 12. 3|) . In § [3l we present our results for each of the 
quasars we analyze. In § [4] we discuss various aspects 
of our results, such as the origin of the constraints we 
obtain, and the robustness of our results in the face of 
uncertainties in the quasar template spectrum. In §(5[ we 
summarize our key findings and present our conclusions. 

2. ANALYSIS 

In this section, we describe the three components of 
our analysis: (i) a hydrodynamical simulation to model 
the absorption (§ 12. 1|) : (ii) modeling the quasars' intrinsic 
emission spectrum (§ 12.2(1 ; and (iii) the statistical com- 
parison of observed and simulated spectra (§ 12. 3p . Unless 
stated otherwise, we adopt the background cosmologi- 
cal parameters (^a, ^m, ^b, n, er 8 , H ) = (0.76, 0.24, 
0.0407, 1, 0.76, 72 km s _1 Mpc" 1 ), consi stent with the 
three- year results by the WMAP satellite (jSpergel et al.1 
[2006ah . We use redshift, z, and the observed wavelength, 
A Q bs, interchangeably throughout this paper as a mea- 
sure of the distance along the line of sight. Unless stated 
otherwise, we quote all quantities in comoving units. 

2.1. Simulating Absorption in the IGM 
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In order to simulate the absorption spectra, we make 
use of a hydrodynamical simulation box in a ACDM uni- 
verse at redshift z = 6. The det ails of the simulation 
are described in ICen et alj (|2003f) . and we only briefly 
discuss the relevant parameters here. The box is 11 hr x 
Mpc on each side, with each pixel being about 25.5 hr x 
kpc. This scale resolves the Jeans length in the smooth, 
partially-ionized 3 IGM by more then a factor of 10. The 
box also has a slightly different set of cosmological pa- 
rameters: (fi A ,^M,^6,n,CT 8 , H ) = (0.71,0.29,0.047, 1, 
0.85 4 , 70 km s^ 1 Mpc^ 1 ). 

Density and velocity information was extracted from 
the simulation box along randomly chosen lines of sight. 
The step size along each LOS was taken to be 5.1 
h^ 1 kpc, which resolves the Lyman a Doppler width 
in the partially-ionized IGM by more than a factor of 
40. The exact value of the step size was chosen some- 
what arbitrarily, and does not influence the results as 
long as it adequately resolves the Doppler width. At 
each step, the density and velocity values were aver- 
aged for the neighboring pixels, and weighted by the 
distance to the center of the pixels. We extended each 
LOS by the common practice of randomly choosing a 
LOS thro ugh the b ox, and stacking sev e ral segments 
together (ICenl Il994t iMesinger et ail 120041 : iKohler et all 
120051 iBolton fc HaehneltJ |2006[ K The r.m.s. mass fluc- 
tuation on the scale the box is a(z — 6.3) ~ 0.15. 
This is a small enough value that neglecting wavelengths 
longer than the box size is justified in our analysis. 
Note however, that this is not the case for modeling 
highly non-linear proc esses, such as reionization (e.g. 
Bar kana fc Loebll2004bD . When using such a LOS in sim- 
ulating a given quasar spectrum, density values from the 
simulation were scaled as (1 + z) 3 . 

The total optical depth due to Lyman a absorption, 
t, between an observer at z = and a source at z = zq, 
at an observed wavelength of A b s = Ao(l + zq), is given 
by: 



T(A ob 



:q dt 

dz c— n{z) xm(z) u 
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where n(z), v(z), xm(z) are the total hydrogen number 
density, the component of the peculiar velocity along the 
LOS, and the local hydrogen neutral fraction, respec- 
tively, at redshift z. The Lya absorption cross section, 
to first order in v(z)/c, is given by er[A bs/(l + z )/(l — 
v(z)/c)}. 

3 The temperature of a neutral IGM before any ionization would 
be several orders of magnitude lower than in the partially-ionized 
case, and the Jeans length in the IGM would be under- resolved by 
a factor of ~6. 

4 The three-y ear WMAP results yielded a smaller value of 
ct 8 = 0.76 ± 0.05 llSpergel et alj|2006 b). A smaller erg would result 
in a narrower distribution of optical depths than those obtained 
from the box we use. Although this does not exactly mimic the 
effect of increased damping wing absorption from a more neutral 
universe, which shifts the optical depth distribution to higher val- 
ues, it is likely that if confirmed, a smaller as would weaken our 
lower limits on the IGM neutral fraction below. Estimates for the 
sensitivity of our results on our choice of cosmological parameters 
is quantitatively difficult without running a suite of models. We 
note however that analysis combining the three-year WMAP re- 
sults with data fr om the Lym an— a forest suggests a higher value of 
erg = 0.86 ±0.03 (Lewis 2006), which is consistent with our choice. 



Since high-redshift quasars are surrounded by their 
own highly ionized HII regions, it is useful to express 
r as a sum of contributions from inside (tr) and outside 
(t£>) the HII region, r — tr + tjj. Defined as such, the 
resonant optical depth, tr, would be given by Eq. (fTJ), 
with the lower limit of integration changed to the red- 
shift corresponding to the edge of the HII region, z(Rs)- 
The damping wing optical depth td, would be given by 
Eq. (QT) , with the upper limit of integration changed to 
z(R s p. 

Since cosmic hydrogen is highly ionized at low red- 
shifts, one can replace the lower limit of integration in 
Eq. (fT]) with z on( j, denoting the redshift below which HI 
absorption becomes insignificant along the LOS to the 
source. In our analysis, we chose the value of z en d sep- 
arately for each quasar to correspond approximately to 
the blue edge of its GP trough. This roughly translates 
to z G nd ~ 0.2-0.3 for our three quasars. We find that the 
exact choice of z en d is not important, because most of the 
contribution to to comes from neutral hydrogen close to 
z(Rg). Inside the HII region, we assume an IGM tem- 
perature of T = 2 x 10 4 K. Outside the HII region, for 
the case of a mostly ionized universe, we assume an IGM 
temperature of T — 10 4 K, while the temperature in a 
neutral universe was taken to be T — 2.73 x 151 ( ^f) 2 K, 
valid for z < 150 (Peebles 1993). Our results are fairly 
insensitive to the exact temperature used. The Doppler 
width of the Lyman-a absorption cross section scales as 
vn oc T 1 / 2 , but the total integrated area under the cross 
section is independent of temperature. 

Assuming ionization equilibrium, we compute the neu- 
tral hydrogen fraction at each point along the LOS with: 



^Hl(rQ + Tbg) = n 2 {l - XYLlfa-B 



(2) 



where Tq and Tbg are the number of ionizations per 
second per hydrogen atom due to the quasar's and the 
background's ionizing fluxes, respectively, and as is the 
case-B hydrogen recombination coefficient. The LHS of 
Eq. ([2]) accounts for the number of hydrogen ionizations 
per cm 3 , while the RHS accounts for the number of hy- 
drogen recombinations per cm 3 . Solving Eq. ([2]) for xhi 
one obtains: 



-b-V& 



(3) 



where b = —2 - (Tq + T B c)/(naB)- We further 
parametrize the quasar's ionization rate with: 



Q 



(4) 



where T Q = (4 7 rr 2 )- 1 /~ 1.55 x 10 31 [v/u H )~ l ' & [(1 + 
z)/(l + z Q )}-°- 8 (o-/hv)dJ's~\ with v H = 3.29 x 10 15 Hz 
being the ionization frequency of hydrogen and r being 
the luminosity distance between the source and redshift 
z. Tq results from redshifting a power-law spectrum 
with a slope of i/~ 1,8 , normalized such that the emission 

5 Note that the terms "resonant" and "damping wing" are only 
appropriate in describing the contributions of tr and To inside 
the HII region (i.e. at A t, s > A a [l + z(R$)]). Outside of the 
HII region (A oba < A a [l + z (/£.<,•)]). tjj actually integrates over the 
resonant part of the absorption cross section and tr integrates over 
the damping wing (see Fig. [TJ. 
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rate of ionizing photons per se cond is 1.3 x 10 57 , ma tch- 
ing the cmissivity one obtains (Haiman fc Cenll200l for 
SPSS J1030+052 4 using a standard template spectrum 
(El vis et al.lll994f ). Note that the normalization infer red 
from the more recent template of iTelfer et al.l (|2002j ) is 
a factor of of ~ 5 higher. This roughly translates to an 
emission rate of ionizing photons of 

N Q = f r x 1.3 x 10 57 s- 1 . (5) 

We ignore radiative transfer along the LOS (and also 
the ionization of helium). Instead, outside of the quasar's 
HII region, we compute xm assuming the quasar con- 
tributes nothing to the ionizing flux, i.e. Tq(R > Rg) = 
0, while inside the HII region, we assume the gas is opti- 
cally thin in the ionizing continuum. This assumption of 
a sharp HII region boundary is an approximation, valid 
in the regime where either the quasar's spectrum is soft, 
or where a strong soft ionizing background flux domi- 
nates over that of the quasar already at radii smaller 
that the mean free path of the typical ionizing quasar 
photon ( as is the case of th e proximity effect at lower 
redshift; iBaitlik et ail 119881) If the quasar has a hard 
spectrum, and the background intensity is low, then the 
edge of the HII region will be "blurred" . A similar blur- 
ring would arise if additional ionizing sources inside the 
HII region, with an extend ed spatial distribution , would 
contribute significant flux ([Wvithe fc Loebll2006l ), a pos- 
sibility we do not address here. However, already strong 
upper limits can be placed on the extent of such blur- 
ring from the offset between the wavelengths of the Lya 
and Ly/3 GP troughs. Since Ly/3 absorption is less sen- 
sitive than Lya, the amount of offset between the GP 
troughs gives a measure of the spatial (i.e. wavelength) 
evolution of r close to the edge of the HII region, with 
a large offset indicating slow r ev olution, and a small 
offset indication rapid r evolution. iMesinger fc Haimanl 
(|2004D studied this for the case of J1030+0524, but even 
stronger constraints could be placed on the spectra of 
J1148+5251 and J1623+3112, since they do not show 
any offsets between the GP troughs (Kramer et al., in 
preparation). Furhtermore, neglecting possible obscura- 
tion effects arising from optically thick gas and dust close 
to the quasar can bias our determination of the quasar's 
luminosity to lower values. This means that our deter- 
mination of the quasar's luminosity should be viewed as 
the observed luminosity along the LOS, and equal to the 
intrinsic luminosity in the low opacity limit. 

Finally, we define the neutral hydrogen fraction in the 
IGM, x^ M , with Eq. ([3]) using the mean density of the 
universe at the edge of the HII region, n — no[l+z(Rs)] 3 , 
and Tq = 0. Although Tbg is the more physically rel- 
evant quantity and we use it in our analysis, henceforth 
we will use x^f M as the proxy for Tbg m order to match 
convention. As defined above, x}§ M represents the neu- 
tral fraction at the mean density. The more common 
representation of the neutral fraction in the literature is 
in the form of the volume (or mass) weighted mean, xhi- 
This definition, of course, depends on the assumed den- 
sity distribution. Hence, for comparison purposes, we 
express our results both in terms of x]^ M and xhi, with 
the latter obtained by averaging over our simulation box 
densities and assuming ionization equilibrium. 

To summarize, our analysis has three free parameters: 



1. Rs, the radius of the quasar's HII region, 

2. xj^f M , the neutral hydrogen fraction in the IGM at 
mean density, 

3. /r, the quasar's ionizing luminosity, in units of 
1.3 x 10 57 s" 1 . 

To obtain an intuition for the effect of varying these free 
parameters, in Figure [T] we plot tr (solid curve) and 
td (dashed curve) for a typical simulated LOS, for pa- 
rameter choices (R s , %m M , fr) = (40 Mpc, 1,1). As 
stated above, the total Lyman a optical depth is the 
sum of these two contributions, r = tr + tjj. Roughly 
speaking, changing Rs moves the dashed (td) curve left 
and right, while changing x^f M moves this curve up and 
down in Figure [TJ Changing fr approximately corre- 
sponds to moving the solid (tr) curve up and down. The 
dotted line demarcates roughly the maximum total op- 
tical depth, r ~ 5.5, at which a signal can be detected 
with a confidence greater than l-cr, in the spectral region 
just blueward of the Lya emission line (see discussion in 
§ 12. 3p . Finally, after the optical depth and the associ- 
ated transmission (e~ r ) is calculated according to the 
procedure outlined above, the transmission is averaged 
over w 3.3 A wavelength bins to match the Keck ESI 
spectra with which the simulated spectra are compared 
(note that the actual smoothing length AA bs differs by 
a few percent for each source, since they are at different 
rcdshifts). 

The free parameter Rs was varied in linear increments 
of ARs — 2 Mpc, with the lower limit set by the red edge 
of the Lya GP through in the observed spectra. The 
other parameters were varied in the ranges: 3.2 x 10 -4 < 
£ffl M < 1, in multiples of 5, and 0.1 < fr < 10, in linear 
increments of Afr = 0.3 between 0.1 < /r < 2 and 
A/ r = 2 between 2 < f T < 10. 

2.2. Modeling the Intrinsic Emission Spectrum 

As observational input in our analysis, we make 
use of Keck ESI spectra of the three highest redshift 
quasars discovered to date: J1148+5251 (zq = 6.42), 
J1030+0524 (z Q = 6.28), and J1623+3112 (z Q = 6.22). 
These quasars are the only quasars discovered to date 
with zq > 6.2, as well as the only quasars (aside from 
the BAL quasar J1048+4637) exhibiting complete GP 
trough s, with no detecta ble flux over a wide wavelength 
range (|Fan et al.ll2006b! ). Readers interested in the de- 
tails of the corresponding observations are encouraged to 
consult iFan et alj (|2006bh and lWhite et al.l (|2003f) . 

In order to compare the quasars' observed spectra with 
our simulated spectra, we must know the intrinsic emis- 
sion spectrum from each of the observ ed sources. Obser- 
vations of lower redshift quasars (e.g. ITelfer et al.1 12002) 
suggest that their intrinsic continuum emission is well fit 
by a broken power-law, with a mean spectral slope of 
a ~ —0.7 (f v oc v a ) red ward of the Lya line, and a slope 
of a ~ —1.8 blueward of the Lya line. The high-redshift 
SDSS quasars do not show any evidenc e for deviating 
from this trend (e.g. I White et al.ll2003l ). The precise 
location of the wavelength of the break in the power- 
law is difficult to determine for any given source, due 
to a strong broad Lya emission line superimposed on 
the continuum spectrum. However, for our purposes, we 
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Fig. 1. — Model from a hydrodynamical simulation for the 
optical depth contributions from within (tr) and from outside (tjj) 
the local ionized region for a typical line of sight towards a zq = 
6.42 quasar. This figure was created with parameter choices (Rs, 
X HI M ' /r) = (40 Mpc, 1, 1). The dashed curve corresponds to tx>, 
and the solid curve corresponds to t^j. The total Lyman a optical 
depth is the sum of these two contributions, r = tr + tx>. The 
dotted line demarcates roughly the maximum total optical depth, 
t ~ 5.5, at which a signal can be detected with a confidence greater 
than 1-ct, in the spect ral re gion just blueward of the Lya emission 
line (see discussion in § 12,31 1, In our analysis, the transmission (e _T ) 
is averaged over ~ 3 - 4 A wavelength bins (1/10 of the resolution 
in the figure above) , thus decreasing the fluctuation in the effective 
optical depth. For reference, the redshifted Lya wavelength is at 
9020 A, off the plot. 

are interested in the optical depth at wavelengths within 
just a few tens of rest-frame Angstroms blueward of the 
Lya line center, where the Lya emission line is still dom- 
inant over the continuum. Therefore, we do not model 
the break, and use a single power-law continuum compo- 
nent instead in our emission template. Throughout our 
analysis, we keep the logarithmic slope a = —0.7 fixed 
(as we will demonstrate below, our results are insensitive 
to this choice). 

We derive our constraints below entirely from the blue 
side of the Lya line in our analys is, where the effects o f 
the IGM can be felt most strongly (jMesinger et alJl2004f ). 
To this end, we make use of "clean" spectral regions on 
the red side of the Lya line, where the resonant attenua- 
tion is negligible (tr ~ 0), and the damping wing atten- 
uation is less than 10% (eT TD > 0.9) even in the case of 
a fully neutral IGM and a small Rs- We fit the NV emis- 
sion line at (1 + zq) 1240.81 A with a single gaussian, and 
the continuum with a single power-law. Additionally, we 
fit the Lya line with a single gaussian for J1148+5251 
(with a fitted width of AA obs = H3A) and J1623+3112 
(with a fitted width of AA Q b s = 8lA) and with a sum of 
two independent gaussians for J1030+0524 (with fitted 
widths of AAobs = 45A and AA D b s = 113A). We have 
found that adding a second gaussian does not improve 



the fits to the Lya line of J1148+5251 and J1623+3112. 
However, the Lya line of J1030+0524 requires two gaus- 
sians for an accurate fit, as it contains a broad and nar- 
row component, similar to the sha pe of a Voigt profile 
(|Ho et alJH997LlGreene fc Hdl200l . The redshifted line 
centers of the gaussians were held fixed in the fitting pro- 
cedure, set to (1 + zq)X , where A = 1215. 67A for Lya 
and A = 1240.81A for N V. The height and width of the 
gaussians were then found by fitting the flux in pixels on 
the red side of the line center. Care was taken, in addi- 
tion, to excise the spectral region immediately redward 
of the Lya line center, as this region could be suscepti- 
ble to non-negligible absorption by dense infalling gas in 
the foreground, close to the q uasar (e.g. iMesinger et al.1 
120041: IBarkana fc Loebll2004al ). Indeed, the quasar spec- 
tra do exhibit sharp drops in the observed flux as one 
approac hes the line center from the red side (see Fig. [2j 
see also IBarkana fc Loebll2003l ). This region was excised 
by eye from the fitting procedure. The amplitude of the 
continuum was found separately, by fitting the flux in re- 
gions of the spectrum further to the red, that are a-priori 
known not to suffer metal line absorption, corresponding 
to roughly 1300 - 1350 A in the rest frame. We explic- 
itly avoid spectral sub-regions with obvious emission or 
absorption lines. More specifically, the spectral regions 
in the observed frame used in the fit for the continuum 
were (demarcated with horizontal lines in Figure [2]): 

• J1148+5251: 9100-9125A, 9128-9190A, 9220- 
9500A, 9700-10000A; 

• J1030+0524: 8880-8990A, 9020-9100A, 9840- 
9930A, 9960-10020A; 

• J1623+3112: 8830-8895A, 8980-9040A, 9790- 
9863A, 9880-9907A; 

In Figure [21 we show the data for the observed spectra, 
F'(Aobs), with 1-ct error bars as well as the continuum + 
Lya + N V fit to the intrinsic emission, Fo(A bs), gen- 
erated by our procedure (solid curve). The Lya optical 
depth can then be inferred in each pixel separately as: 
7"obs(A bs) = -ln[F(A bs)/^u(Aobs)]- 

2.3. Comparing Simulated and Observed Spectra 

After the observed optical depths, T b s (A b s ), are ob- 
tained as discussed in § 12.21 and the simulated optical 
depths, r s i m (A bs), are obtained for every LOS and point 
in our parameter space as discussed in § 12.11 we wish 
to statistically compare the simulated and observed op- 
tical depths for each quasar. More precisely, for each 
quasar and every point in our three-dimensional param- 
eter space, (Rs, i™! /r)j we want to answer the ques- 
tion: what is the likelihood that the list o/r b s (A bs) val- 
ues and the list o/r S i m (A b s ) values were drawn from the 
same underlying probability distribution? Note that we 
have a fixed list of about 30 r bs(A bs) values for each 
quasar, but we generate a different list of 3000 T S i, n (A b s ) 
values, combining 100 LOSs. for each choice of (Rs, 

4! M , fr). 

We use the K olmogorov-Smirnov (K-S) test (e.g. 
iPress et alJ H992) to compare the cumulative probability 
distributions (CPDFs) - i.e., histograms - of r b s (A b s ) 
and T s i m (A b s )- The usefulness of the K-S test is that 



() 




Fig. 2. — Observed spectra (shown as points with 1-cr error bars) and the continuum + Lya + N V fit of the intrinsic emission spectrum 
generated by our procedure (solid curve). The spectral regions (on the red side of the Lyo line center) that we used in the fit are marked 
with horizontal lines. The inlaid panels zoom in on the spectral region (on the blue side of the Lya line center) that were then used to 
compare the absorption optical dept h r with hydrodynamical simulations. The vertical dashed lines bracket the wavelength ranges used 
to construct r histograms (see The spectra shown are those of J1148+5251, J1030+0524, J1623+3112 (left to right). 



it provides a figure of merit which can be directly in- 
terpreted as the likelihood that two distributions were 
drawn from the same underlying distribution, without 
the need for the ad-hoc binning of data, and the resulting 
loss of information, required by such methods as the % 2 
test. The drawback of using the K-S test is that comput- 
ing absolute confidence levels around our best fit values 
would require prohibitively expensive Monte-Carlo sim- 
ulations. We caution the reader that, as a result, the 
relative likelihoods we quote below do not represent tra- 
ditional "error-bars" . 

Our sta tistical analysi s is qu ite similar to that cm- 
ployed by iRollinde et al.l (|2005[ ) to model density struc- 
ture surrounding moderate redshift (z ~ 2) quasars. The 
main difference in technique is that wh ile our present 
appro ach uses a parameterized model, IRollinde et al.1 
(2005) fit for the intrinsic spectrum directly from the 
low-r regions of the spectrum. Using simulated spec- 
tra, they show that by applying the K-S test on optical 
depth CPDFs, one can statistically discriminate among 
different models of the di s tance -evolution of the optical 
depth. In iMesinger et al.l (|2004D . we use a similar statis- 
tical test to show that the moderate degeneracy between 
Rs and can be broken with a single spectrum, if 

Rs is moderately constrained. This prior analysis did not 
take into account the additional constraint arising from 
the evolution of the CPDF with distance, as our current 
analysis does. Nevertheless, we caution the reader that 
these tests only confirm the ability to statistically recover 
input parameters using our statistical analysis. They do 
not test the validity of our model. 

The optical depth CPDF should have a strong depen- 
dence on wavelength (since the gas is more highly ion- 
ized near the quasar). In principle, the information con- 
tained in this dependence could be extracted and used to 
further constrain model parameters. For example, one 
could compare the CPDFs of T b s (A b s ) and r S i m (A b s ) 
in several wavelength bins for each quasar. Unfortu- 
nately, in practice, one also needs to have many points 
in each individual bin, to be able to sample the under- 
lying distribution in that bin. In one extreme, as the 
wavelength bin size approaches the pixel width, the in- 



formation from the wavelength dependence is fully pre- 
served; however, only a single observed pixel is available 
to compare with a corresponding pixel from each LOS, 
rendering the use of the K-S test unreliable. The con- 
verse is true in the other extreme, where the wavelength 
bin size approaches the entire spectral region we analyze: 
the shape of the CPDF is accurately captured, but in- 
formation from the wavelength-dependence is destroyed. 
Here we (somewhat arbitrarily) divide the spectral re- 
gion of analysis into wavelength bins of width ~ 25-30 
A, so that each bin contains between 6-8 pixels. The 
two blue-most bins were chosen bracket the edge of the 
GP trough in the observed spectrum. Furthermore, we 
remove the region 20-30 A blueward of the line center 
from our analysis, corresponding to the expected mean 
radius of the large-scale overdense re gion surrounding 
such quasars (Ba rkana fc Loebl l2004af) . This step was 
necessary, as our simulation box size is too small to con- 
tain a statistical sample (or even a single occurrence) of 
such large overdensities, required for their accurate mod- 
eling. The wavelength bins we use in our analysis are 
shown by the vertical dashed lines in the inlaid panels of 
Figure [2] We verify that our results are fairly insensitive 
to our choice of bin size, i.e. the relative K-S likelihoods 
between models change little with decreasing bin size. 

In order to incorporate an uncertainty into our emis- 
sion template, we add Gaussian-distributed, uncorre- 
lated flux variations to each pixel in the template. We 
do this by regarding our fitted value of the template, 
Fo(Xobs), as the mean value of the true flux distribution, 
and the "true" emission at A b s is drawn from a Gaussian 
distribution around that mean, with a standard deviation 
o-ftux(Kbs)- Hence, we replace every T sim (A obs ) value by 
a set of 1000 values, defined by 

r s ? m mplcd (A obs ) = - lnle—^ x A] , (6) 

where A represents our emission uncertainty, and is 
drawn from a Gaussian distribution with a mean of 1, 
and a standard deviation of o~fi ux 

(AobsV-FMAobs). Wc 
set <7/7„ x (A b s )/Fo(A b s ) = °- 3 > wn i cn corresponds ap- 
proximately to the typical quasar-to-quasar scatter in 
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the continuum level imm ediately redward of Lyman a 
(|Vanden Berk et all l200lh . However, we have verified 
that our results are not sensitive to this choice. For ex- 
ample, in the case of J1623+3112, we find that in the 
range 0.1 < o~fi ux (A bs)/-Fo(Aobs) < 0.5, the peak K-S 
probability remains unchanged, and the confidence lim- 
its quoted below change at most by a single pixel in our 
parameter space. 

Finally, in comparing spectra, one needs to define what 
one considers to be "zero flux" , as simulations have near- 
infinite precision and observations do not. To this end, 
we define a maximum allowed value of the optical depth, 
MAX[t] = 5.5. This approximately corresponds to a \-u 
detection in our observed spectra at wavelengths around 
the edge of the GP trough, where the majority of pixels 
with fluxes less than this value are located (see Fig. [2|). 
In our comparison, all values of r > 5.5 are treated the 
same, regardless of their actual value. 

To summarize, for each quasar we divide the spectral 
region roughly blueward of ~ (1 + zq)Aq,-25A and red- 
ward of ~ [1 + z(i?s)]A a -25A into 3-5 bins of approx- 
imately equal wavelength width, containing 6-8 pixels 
each. In each bin we compare the CPDFs of r bs(A bs) 

and T s s ™ ple (A bs) (the latter concatenated over 100 simu- 
lated LOSs). For every point in our 3D parameter space, 
these two distributions are compared via the K-S test, 
which produces a probability to be interpreted as the 
likelihood that the two distributions were drawn from 
the same underlying distribution. Finally, the K-S prob- 
abilities for each wavelength bin are multiplied to give a 
likelihood estimate for that point in the 3D parameter 
space. 

3. RESULTS 

Our results for the likelihoods of various parameter- 
combinations are shown in Figure [3] for J1148+5251 
{top), J1030+0524 (middle), and J1623+3112 (bottom). 
The cross in each case marks the point in the 3D param- 
eter space where the likelihood peaks. In order from 
lighter to darker, the yellow, green, and blue squares 
correspond to points in our parameter space with like- 
lihoods that are a factor of 27, 9, and 3 below the peak 
likelihood, respectively. 6 Parameter combinations with 
likelihoods smaller than l/27th of the peak are shown as 
empty squares. 

Below, we will discuss each quasar in turn. However, 
it is encouraging to note first that the isocontours in 
Figure [3] behave as expected, following the moderate de- 
generacies in the parameters. Namely, since a stronger 
damping wing (td) contribution can be obtained either 
by increasing the neutral fraction or by decreasing the 
size of the HII region, the isocontours should be oriented 
from bottom-left to upper-right in each panel. Similarly, 
since a smaller ionizing flux (i.e. greater tr) can roughly 
be compensated for in the total optical depth, tr + td, 
by increasing the size of the HII region or decreasing 
the neutral fraction (both decrease Trj), the isocontours 
should shift monotonically from the upper-left to the 

6 In order to obtain absolute confidence limits, we would need 
to perform a set of computationally expensive Monte-Carlo simu- 
lations. For simplicity, we avoid such computations, and instead 
contend ourselves by quoting results in terms of the relative offset 
from the peak likelihood. 



bottom-right with increasing fr- Both of these trends 
are evident in Figure[3] For example, the colored squares 
(corresponding to K-S likelihoods that arc within a factor 
of 27 of the peak likelihood) in the middle row cover most 
of the displayed parameter space in the left most panel, 
though there are no pixels with high likelihoods. In- 
creasing the quasar's ionizing flux (i.e. going through the 
panels from left to right), shifts the colored squares in- 
creasingly towards the bottom-right corner of the panel. 

3.1. J1U8+5251 

Likelihood estimates for J1148+5251 (zq = 6.42) can 
be seen in the top row of three panels in Fig. [3l The left, 
center, and right panels correspond to values of /r= 0.7, 
1.6, and 1.9, respectively. The peak likelihood of 3.0% 
occurs at (R s , a;Jg M , /r) = (40 Mpc, 0.2, 1.6). This is the 
smallest peak likelihood we obtain (c.f. peak likelihoods 
of ^30% for the two quasars below), and most likely 
means that our model isn't as accurate in describing the 
physical spectra of this particular quasar. The volume 
averaged neutral fraction at the peak is xhi = 0.16. The 
range of parameter values whose likelihood estimates are 
within a factor of 3 of the peak likelihood is encompassed 
by: 

• 40 Mpc < R s < 42 Mpc, 

• a ; Hi M 1-0 or x m < 1.0, 

• °- 7 < 1-9- 

Note that here as well as below, the end points of the 
quoted ranges correspond to the minimum or maximum 
value of the parameter for which there exists a combina- 
tion of the other two parameters giving a K-S probability 
within a factor of three from the peak probability. Note 
that these ranges are similar in spirit to "marginalized 
errors" , in that the other parameters are allowed to vary, 
rather than fixed (but they are not actual marginalized 
errors; see caveat mentioned in the footnote above). 

3.2. J1030+0524 

Likelihood estimates for J1030+0524 (zq = 6.28) are 
shown in the middle row of panels in Figure [3l The left, 
center, and right panels correspond to values of /r= 0.4, 
1.0, and 1.3, respectively. The peak likelihood of 34% 
occurs at (R s , x^ M , f r ) = (41 Mpc, 1, 1.0). The range 
of parameter values whose likelihood estimates are within 
a factor of 3 of the peak likelihood is encompassed by: 

• 37 Mpc < R s < 45 Mpc, 

• x m M ^ !-° or j hi < 1-0, 

• 0.7 < f r < 1.9. 

As mentioned earlier, J1030+0524 has a unique (at 
least among our sample of three quasars) feature. 
Namely, the Ly/3 GP through is noticeably offset from the 
Lya GP through. The mere presence of flux (irrespective 
of its properties) in the Ly/3 region of the spectrum can 
be used to set a minimum size of the surrounding HII 
region, Rs > 40 Mpc. The region excluded by this prior 
is shaded over in the middle row in Figure [3l Taking this 
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Fig. 3. — Likelihoods for J1148+5251 (top), J1030+0524 (middle), and J1623+3112 (bottom) in the 3-dimensional parameter space of 
(Rs, ^hi M i /r)- Each panel shows contours in a 2D slice of the 3D likelihood surface. The cross marks the parameter combination with 
the peak likelihood value (see § 12.31 for how the likelihoods were derived). In order from lighter to darker, the yellow, green, and blue 
squares correspond to points in our parameter space with likelihoods that are a factor of 27, 9, and 3 below the peak value, respectively. 
Parameter combinations with likelihoods below l/27th of the peak are shown as empty squares. J1148+5251: The peak likelihood of 3.0% 
occurs at (Rs, x hi M i fr) = (40 Mpc, 0.2, 1.6). The left, center, and right panels correspond to values of /r= 0.7, 1.6, 1.9, respectively. 
J1030+0524: The peak likelihood value of 34% occurs at (Rs, z^f M , fr) = (41 Mpc, 1, 1.0). The left, center, and right panels correspond 
to values of /r= 0.4, 1.0, 1.3, respectively. The shaded region shows values ruled out by the presence of flux in the Lyman j3 region, which 
sets the additional constraint R s > 40 Mpc. J1623+3112: The peak likelihood value of 39% is at (R s , z^f M , / r ) = (29 Mpc, 1, 0.7). The 
left, center, and right panels correspond to values of /r= 0.4, 0.7, 1.0, respectively. 



lower limit on Rs into account, the range of values whose 
likelihood estimates are within a factor of 3 of the peak 
likelihood shrinks to: 

• 41 Mpc < R s < 45 Mpc, 

• 0.04 <x^ M or 0.033 <x m , 

• °- 7 %fr < 1-3. 



3.3. J1623+3112 

Likelihood estimates for J1623+3112 (zq = 6.22) are 
shown in the bottom row of panels in Figure [3] The 
left, center, and right panels correspond to values of fr= 
0.4, 0.7, and 1.0, respectively. The peak likelihood of 



39% occurs at (Rs, 



„IGM 



fr) = (29 Mpc, 1, 0.7). The 



range of parameter values whose likelihood estimates are 
within a factor of 3 of the peak likelihood is encompassed 
by: 
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J1623+3112 




Fig. 4. — Distributions of simulated optical depths for the 
three wavelength bins used in our analysis of J1623+3112. The 
solid blue, long-dashed red, short-dashed green, and dotted purple 
curves correspond to parameter combinations of (Rs, X ^ M , fr) = 
(29 Mpc, 1.0, 0.7), (29 Mpc, 0.04, 0.7), (33 Mpc, 1.0, 0.7), and (29 
Mpc, 1.0, 0.1), respectively. The crosses are located at the set of 
observed optical depth values. The solid blue curve corresponds to 
the parameter combination with the peak likelihood value, while 
the parameter combinations represented by the other curves have 
likelihood values below l/27th of the peak. The sharp jump at the 
largest r bin and the overlaying of seven crosses in the top panel are 
due to the finite detection threshold, whereby all values of r > 5.5 
are indistinguishable from one another, and hence fall in the same 
bin in this figure. The figure demonstrates that our lower limit 
on £j£f M comes from the low— r tail predicted in models with low 
neutral fractions or those corresponding to regions of parameter 
space which could mimic a low neutral fraction (see the discussion 
on isocontours in §13} . No points corresponding to these tails are 
seen in the data. Note also that while some models might perform 
well in one wavelength bin, only the solid blue curve is consistent 
with the observed data in all wavelength bins. 



• 25 Mpc < R s < 29 Mpc, 

• 0.04 <x^ M or 0.033 < x m , 

• 0.4 <fr< 1-3. 

4. DISCUSSION 

In this section, we discuss a few important aspects of 
the above results. First, it is encouraging that the sim- 
ple methodology we outlined can deliver statistically in- 
teresting constraints, even from the few dozen available 
spectral pixels. Clearly, it will be interesting to apply 
a similar analysis to a larger sample of sources in the 
future. 

Perhaps our most interesting result is the lower limit 
we obtain, for two of the quasars, on the neutral frac- 
tion. We would like to therefore understand where these 
constraints actually come from. To this end, in Figure 2] 
we show the simulated optical depth PDFs for the three 
wavelength bins used in our analysis of J1623+3112. The 



solid blue, long-dashed red, short-dashed green, and dot- 
ted purple curves correspond to parameter combinations 
of (R s , x^ M , fr) = (29 Mpc, 1.0, 0.7), (29 Mpc, 0.04, 
0.7), (33 Mpc, 1.0, 0.7), and (29 Mpc, 1.0, 0.1), respec- 
tively. The crosses are located at the set of observed 
optical depth values. The solid blue curve corresponds 
to the parameter combination with the peak likelihood 
value, while the parameter combinations represented by 
the other curves have likelihood values below l/27th of 
the peak. The sharp jump corresponding at the largest 
t bin and the overlaying of seven crosses in the top panel 
are due to the finite detection threshold, whereby all val- 
ues of r > 5.5 are indistinguishable from one another, 
and hence fall in the same bin in this figure. The fig- 
ure clearly demonstrates that our lower limit on xj^f 
arises from the low-r tail predicted in models with low 
neutral fractions or those corresponding to regions of pa- 
rameter space which could mimic a low neutral fraction 
(see the discussion on isocontours in § [3]) . While the solid 
blue curves are consistent with the distribution of points, 
these other curves all have tails that are "too long" : no 
points corresponding to these tails are seen in the data. 
Note also that while some models might perform well in 
one wavelength bin, only the solid blue curve is consistent 
with the observed data in all wavelength bins. 

4.1. Sanity Checks 

As our analysis technique is as yet untested, it is use- 
ful to take a step back, and verify that our results make 
sense. As already mentioned, the likelihood iso-contours 
behave as expected. Here we discuss several other en- 
couraging checks on the validity of our analysis. We cau- 
tion the reader, however, that these are not meant to 
be a proof of the accuracy of our analysis, and instead 
should be regarded merely as consistency checks. 

First, the procedure here recovers fairly well the re- 
sult of our previous analys is in the case of J1030+0524 
(jMesinger fc Haimanll2004h . This need not be the case, 
since the nature of the two analysis procedures are dif- 
ferent; namely, our previous work focused only on gross, 
approximate Lya and Ly/3 transmission features near the 
edge of the HII region. Nevertheless, in both cases, the 
peak likelihood occurs at the same value of x^ M and 
Rs, with only the value of fr be ing somewhat larger 
here (/r=1.0 instead of /r=0.4 in iMesinger fc Haimanl 
2004). We attribute this shift in fr primarily to the in- 
clusion of wavelength bins further inward from the HII 
region edge, where the effects of Rs and x^ M are felt 
less strongly and fr wields the most statistical weight. 
In particular, the likelihoods in the two blue-most wave- 
length bins are comparable for fr = 1.0 and fr = 0.4 
(given x^ M — 1.0 and Rs = 41Mpc); however, the like- 
lihoods in the two red-most wavelength bins are ~ 7-8 
times greater for fr = 1.0 than fr = 0.4. 

Second, our results on the neutral fraction agree in 
relative terms with estimates obtain ed from transmis- 
sion in GP troughs (|Fan et al.ll2006al ). J1030+0524 and 
J1623+3112 both have very dark GP troughs in Lya 
and Ly/3, and thus are only able to provide lower limits 
on the optical depth and corresponding neutral fraction; 
however, J1148+5251 has transmission gaps in the GP 
troughs, thus favoring a smaller optical depth and corre- 
sponding neutral fraction. This is qualitatively similar to 
our results: we were unable to constrain x^ M from the 
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spectra of J1148+5251 obtaining only a shallow peak in 
the likelihood at x^ M = 0.2; however, both J1030+0524 
and J1623+3112 had peak likelihoods at x^f M = 1, and 
strongly prefer Xjjf M > 0.04. 

Third, the maximum likelihoods of fr (i.e. the rela- 
tive strength of the quasars' ionizing luminosity) match 
the relative strengths of the quasars' continua redward of 
Lya. The values of fr corresponding to the peak likeli- 
hoods are 1.6, 1.0, and 0.7 for J1148+5251 J1030+0524 
and J1623+3112 respectively. Looking at Figure [2l 
one can note that this is also the order of the relative 
strengths of the quasars' continua redward of Lya, as one 
would expect if the UV power-law indices and amount of 
obscuration didn't vary greatly among the quasars. 

Finally, the quasar with the smallest HII region, 
J1623+3112 with a peak likelihood at R s = 29 Mpc, 
also has the highest favored value of x[jf M = 1.0 and 
the lowest favored value of fr- This makes sense, since 
a fainter quasar, embedded in a more neutral medium 
should have a smaller HII region, given a similar lifetime. 

4.2. Uncertainties in the Intrinsic Emission Spectra 

Although we have added gaussian noise to our inferred 
template spectrum, as described above, we wish to get 
a sense of how correlated errors, modifying the broader 
shape of the template spectrum, would affect our results. 
The general concern is that mis-estimates of the flux, 
correlated over many pixels, could potentially mimic the 
effect of shifts in our parameter space. To get a sense of 
the magnitude of this effect, we chose to vary the width 
of the fitted Lya gaussian for quasar J1623+3112 and see 
how this modifies our results. Naively, one would expect 
a narrower Lya emission line and an accompanying de- 
crease in effective r b s to cause a shift in the parameter 
space to regions favoring smaller optical depths (larger 
Rs, smaller xj^f M , larger fr), with the opposite trends 
for a wider Lya emission line. 

We find that such degeneracies are very weak. For a 
10% narrower line, we find that the peak likelihood re- 
mains at the same parameter-combinations. The range 
of parameter values whose likelihood estimates are within 
a factor of 3 of the peak likelihood changes by at most one 
pixel in our 3D parameter grid: 25 Mpc < Rs < 29 Mpc; 
0.008 < a;^f M < 1.0; 0.4 < fr < 1.6. This is a rather 
small change, and we emphasize that a 10% narrower 
Lya line would actually be an unacceptably poor fit to 
the (red side of the) observed spectrum. A wider Lya line 
is more physically motivated, since strong damping wing 
absorption might be able to somewhat attenuate flux far 
into the red side of the line. Hence, we vary the width 
of the Lya line, well in excess of the absorbing potential 
of a strong damping wing. Again, we find little change 
in our results. Even for a 50% wider line, the param- 
eter combination at peak likelihood remains the same, 
while the range of parameter values whose likelihood es- 
timates are within a factor of 3 of the peak likelihood 
changes by at most one pixel to: 25 Mpc < Rs < 31 
Mpc; 0.2 < x*g M < 1.0; 0.4 < fr < 1.3. As expected, an 
intrinsically wider Lya line shows a slight preference for 
a more neutral universe. 

Overall, the above exercise suggests that our results 
are robust and conservative, at least to mis-estimates of 



the Lya line width. This does not, of course, preclude 
the possibility that the intrinsic Lya lines of the z w 6 
quasars possess some deviations from a gaussian shape. 
For example, if there is an asymmetry, in the sense of 
a steeper drop of flux on the blue side of the line than 
on the red side, our technique would have mistakenly 
attributed this steeper drop to absorption. This implies 
that we might have overestimated the neutral fractions. 
At present, there is, however, no evidence of such an 
asymmetry in the emission lines of low redshift quasars 
(where both sides of the lines can be seen unabsorbed) , 
nor is there evidence, from the red side of the lines, that 
the high-z quasars' lines have shapes different from those 
of their lowcr-redshift counterparts. 

Finally, we also note that our spectral fits, describing 
the intrinsic emission line on the unabsorbed red side 
of the line, leave typical flux residuals of only ~ 1%. 
Similar residuals are presumably present on the blue side 
of the line, as well. However, the absorption exp(— td) 
caused by the IGM damping wing on the blue side of the 

line (see Fig. [T]), has an amplitude of ~ (0.1 1)xhi, 

exceeding the residuals for xhi > 0.01. Furthermore, in 
order for the residuals to mimic the IGM damping wing, 
they would have to be coherent (monotonically increasing 
with wavelength). 

4.3. Biases in Local Density and Velocity Fields 

As mentioned above, we remove the region 20-30 A 
blueward of the line center from our analysis, correspond- 
ing to the expected mean radius of the large-scale over- 
dense region surrounding such quasars (|Barkana fc Loebl 
l2004al) . This is a necessary step as present-day simu- 
lations cannot statistically model such rare, large-scale 
overdensities. However, if the density bias extended 
into our region of analysis, our inferred x}§ M lower lim- 
its would be conservative and our inferred fr values 
would most likely be underestimated. If the source is 
hosted by a large-scale overdensity, then density would 
decrease with decreasing A b s (increasing distance from 
the source). If such a decreasing density extended into 
our region of analysis, it would translate into a shallower 
rise in tr oc n 2 /T. In matching the total optical depth 
T hya = tr + td, our analysis, which assumes constant 
mean density, would translate this shallower rise in tr 
to a shallower rise in tjj, i.e. a lower value of x^ M than 
appropriate if the density bias was taken into account 7 . 
Note that our results on x[jf M are driven by the wave- 
length dependence of r; a constant offset in the density 
field can be absorbed by a shift in Tq (i.e. fr), which 
dominates over Tbg throughout most of our region of 
interest and is predominantly set by the amount of ab- 
sorption in bins further from the HII region edge where it 
has the most statistical weight. Furthermore, it is possi- 
ble that an infalling IGM could introduce a similar bias, 
but the effect is likely to be small for a bright quasar and 
far away from the line center on the blue side. More- 
over, the effect from the infall would have the same sign 
as the density bias above, making the flux fall off on 

7 Note that a shallower rise in Tp could also be achieved if we 
were overestimating Rg, instead of underestimating a;j^| M . How- 
ever, our inferred Rg values are already at the lowest limit set by 
the onset of GP troughs for two of our quasars and thus can not 
be decreased further. 
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the blue side of the lin e less steeply (|Diikstra et al.ll2006l ; 
iBarkana fc Loe b 2003)). In other words, by ignoring the 
possible density and velocity biases, we might only be 
slightly underestimating the neutral fraction. 

4.4. Uncertainties in Intrinsic Redshift 

A different source of possible error, in estimates of Rs, 
can result from an incorrect redshift determination. The 
rcdshifts we quote above were determined from high- 
ionization lines such as C IV, Si IV and N V. However, 
high-ionization emission lines can have non-negligible off- 
sets from low-ionization, narrow, or molecular lines in the 
host galaxy. Recently, redshift determinations based on 
such lines have been made for the three qu asars used in 
our a nalysis: z C o = 6.42 for J 1148+525 1 (IWalter et al. 
2003t); ZMgii = 6.31 for J1030+0524 (jlwamuro et al. 
2004); z Mg ii = 6.22 for J1623+3112 (L. Jiang et al. 2006, 
in preparation). These redshift determinations agree 
with the ones based on high-ionization lines used in this 
work for J1148+5251 and J1623+3112; however the Mg 
II redshift is 0.03 larger than the high-ionization line red- 
shift for J1030+0524. If one takes the Mg II redshift to be 
the true systematic redshift, this would roughly increase 
all R s determinations for J1030+0524 by - 30% The 
majority of the effect would be confined to the increase 
in Rs, since the systematic redshift doesn't impact the 
Lya emission fitting, and the isocontours (as seen in the 
middle panels of Fig. [3|) are steep (i.e. a large increase in 
Rs can be compensated for by a small increase in £^f M , 
especially for the lower confidence contours). Since in- 
creasing Rs can roughly be compensated for by increas- 
ing XjjJ or decreasing /p, the secondary effects of such 
a redshift underestimate should be a shift in the favored 
parameter values towards larger x^ M and/or smaller /p. 

5. CONCLUSIONS 

We have made use of hydrodynamical simulations of 
the IGM to create model quasar absorption spectra, and 
compared these with observed Keck ESI spectra of three 
z > 6.2 quasars: J1148+5251 (z = 6.42), J1030+0524 
(z=6.28), and J1623+3112 (z=6.22). We compared the 
CPDFs of observed Lya optical depths to those gener- 
ated from the simulation by assuming various values for 
the size of the quasar's surrounding HII region, Rs (in 
comoving Mpc) , volume- weighted neutral hydrogen frac- 
tion in the ambient IGM, Xm, and the quasar's emission 
rate of ionizing photons, Nq (in 10 57 s -1 ). This ap- 
proach has not yet been applied to the spectra of the 
z > 6 quasars, and we have emphasized some impor- 
tant caveats, most importantly our neglect of radiative 
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transfer effects, and our assumption that the intrinsic 
emission line of the quasars is symmetric around the line 
center. Since large sightline-to-sightline fluctuations in 
IGM properties are likely at these high-redshifts, we an- 
alyzed each quasar independently. We found that our 
approach can constrain Rs to within ~ 10%, and Nq to 
within a factor of ~ 2 for all three quasars. Fairly strong 
constraints on xm are obtained from J1030+0524 and 
J1623+3112, while xm is unconstrained by J1148+5251. 
Specifically, our results are as follows. 

,11148+5251: the parameter combination with the 
maximum likelihood is (Rs, xm, Nq) = (40, 0.16, 2.1); 
the range of parameter values whose likelihood estimates 
are within a factor of three of the peak is encompassed 
by: 40 < R s < 42, x m < 1, 0.9 < N Q < 2.5. 

J1030+ 0524'- the parameter combination with the 
maximum likelihood is (Rs, Shi, Nq) — (41, 1, 1.3); 
the range of parameter values within a factor of three of 
the peak likelihood is encompassed by: 41 < Rs < 45, 

0.033 < x m , 0.9 <N Q < 1.7. 

J 1623+3 112: the parameter combination with the 
maximum likelihood is (Rs, xm, Nq) = (29, 1, 0.9); 
the range of parameter values within a factor of three of 
the peak likelihood is encompassed by: 25 < Rs < 29, 

0.033 < x m , 0.5 <N Q < 1.7. 

It is encouraging that the simple methodology we have 
utilized can deliver statistically interesting constraints, 
even from the few dozen available spectral pixels. Our 
method is different from previous analyses of the GP ab- 
sorption spectra of these quasars, and the lower limits 
we found on the neutral fraction of the IGM surround- 
ing two of the quasars strengthen the evidence that the 
rapid end-stage of reionization is occurring near z ~ 6. 
Clearly, it will be interesting to apply a similar analysis 
to a larger sample of sources in the future. 
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